      subroutine Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid, &
         jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, sdlt, strait&
         ,dx,dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y&
         , c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, vvol, &
         vvolaxis, q, qtilde, aqtilde, ht, itmax, error, srce, rdt, scalsq, wk&
         , phibc, ex, ey, ez, dive, residu, aq, phi, diag, wkl, wkr, wkf, wke, &
         periodicx)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      use boundary

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ncells
      integer  :: nvtxkm
      integer  :: nvtx
      integer  :: itdim
      integer  :: iwid
      integer  :: jwid
      integer  :: kwid
      integer  :: ibp1
      integer  :: jbp1
      integer  :: kbp1
      integer , intent(in) :: itmax
      real(double)  :: tiny
      real(double) , intent(in) :: relax
      real(double)  :: cdlt
      real(double)  :: sdlt
      real(double)  :: strait
      real(double)  :: dx,dy,dz
      real(double) , intent(in) :: error
      real(double)  :: rdt
      real(double)  :: scalsq
      real(double)  :: wk
      real(double)  :: phibc
      real(double)  :: wkl
      real(double)  :: wkr
      real(double)  :: wkf
      real(double)  :: wke
      logical , intent(in) :: precon
      logical  :: periodicx
      integer  :: ijkcell(*)
      integer  :: ijkvtx(*)
      real(double)  :: c1x(*)
      real(double)  :: c2x(*)
      real(double)  :: c3x(*)
      real(double)  :: c4x(*)
      real(double)  :: c5x(*)
      real(double)  :: c6x(*)
      real(double)  :: c7x(*)
      real(double)  :: c8x(*)
      real(double)  :: c1y(*)
      real(double)  :: c2y(*)
      real(double)  :: c3y(*)
      real(double)  :: c4y(*)
      real(double)  :: c5y(*)
      real(double)  :: c6y(*)
      real(double)  :: c7y(*)
      real(double)  :: c8y(*)
      real(double)  :: c1z(*)
      real(double)  :: c2z(*)
      real(double)  :: c3z(*)
      real(double)  :: c4z(*)
      real(double)  :: c5z(*)
      real(double)  :: c6z(*)
      real(double)  :: c7z(*)
      real(double)  :: c8z(*)
      real(double)  :: vol(*)
      real(double)  :: vvol(*)
      real(double)  :: vvolaxis(*)
      real(double) , intent(inout) :: q(*)
      real(double)  :: qtilde(*)
      real(double)  :: aqtilde(*)
      real(double)  :: ht(*)
      real(double)  :: srce(*)
      real(double)  :: ex(*)
      real(double)  :: ey(*)
      real(double)  :: ez(*)
      real(double)  :: dive(*)
      real(double)  :: residu(*)
      real(double) , intent(inout) :: aq(*)
      real(double)  :: phi(*)
      real(double)  :: diag(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: nstart, n, k, j, i, ijk, iter, ijkmax
      real(double) :: lambda, rsum, rsqsum_old, rsqsum_new, bottom, rnorm, &
         dxnorm, xnorm, alfa,  dumdt, dumscal, phibcold
!-----------------------------------------------
!
!
!
       q(ijkcell(:ncells))=0.0d0
       aq(ijkcell(:ncells))=0.0d0

       qtilde(ijkcell(:ncells))=0.0d0
!
      phibcold=phibc
      call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkl, wkr, wkf, wke, phibc&
         , phi, periodicx)
!
!     CALCULATE THE INITIAL RESIDUAL ERROR
!
      call ResidueCG (ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibp1, jbp1&
         , tiny, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, &
         c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, phi, &
         ex, ey, ez, vvol, dive, srce, rdt, scalsq, residu)
!
!       phi(:itdim)=0.0d0
!
      rsum=0.0d0
      rsqsum_new=0.0d0
      do n=1,ncells
         rsum = rsum+ abs(residu(ijkcell(n)))
         rsqsum_new=rsqsum_new+residu(ijkcell(n))*residu(ijkcell(n))
      enddo
!
!      write(*,*)'PoissonCG:  dnorm before=',rsqsum_new
      iter=0
      if (rsum < error) go to 2
!
!     ****************************************************************
!
!     begin conjugate gradient iteration
!
!     ****************************************************************
!
      do iter = 1, itmax
!
      if(iter < 2) then
         do n=1,ncells
            q(ijkcell(n))=residu(ijkcell(n))
         enddo
      else
         rsqsum_old=rsqsum_new
         rsqsum_new=0.0d0
         do n=1,ncells
            rsqsum_new=rsqsum_new  &
              +residu(ijkcell(n))*residu(ijkcell(n))
         enddo
         lambda=rsqsum_new/rsqsum_old
         do n=1,ncells
            q(ijkcell(n))=residu(ijkcell(n))  &
               +lambda*q(ijkcell(n))
         enddo
      endif
!
      phibc=0d0
      call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkl, wkr, wkf, wke, phibc&
         , q, periodicx)
!      write(*,*) 'PoissonCG: iter, rsqsum_old,rsqsum_new=',iter,rsqsum_old,rsqsum_new
!
!     calculate residue
!
!         dumdt=0.0d0
!         dumscal=-1.d0
         dumdt=rdt
         dumscal=scalsq
!
         call ResidueCG (ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibp1, &
            jbp1, tiny, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, &
            c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, &
            vol, q, ex, ey, ez, vvol, dive, qtilde, dumdt, dumscal, aq)
!
       bottom=0.0d0
       do n=1,ncells
          bottom=bottom+q(ijkcell(n))*aq(ijkcell(n))
       enddo
       alfa=rsqsum_new/(bottom+1.d-20)
!
!
!

!
        rnorm=0.0d0
        dxnorm=0.0d0
        xnorm=0.0d0
        do n=1,ncells
           phi(ijkcell(n)) = phi(ijkcell(n))   &
            + alfa*q(ijkcell(n))
           residu(ijkcell(n)) = residu(ijkcell(n))   &
            - alfa*aq(ijkcell(n))
!
           rnorm = rnorm+abs(residu(ijkcell(n)))
           dxnorm = dxnorm+abs(q(ijkcell(n)))
           xnorm = xnorm+abs(phi(ijkcell(n)))
        enddo
!
         if (dxnorm<=error*xnorm .and. rnorm<=error*rsum) go to 2
!
      end do
!     iteration failed
      write (*, *) 'PoissonCG:  iteration relies on the good luck'
!
      return
!
    2 continue
      write (*, *) 'PoissonCG:  iteration converged, iter=', iter
!
!     correct error in the divergence of e0
!
      phibc=phibcold
      call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkl, wkr, wkf, wke, phibc&
         , phi, periodicx)
!
!
!     CALCULATE THE FINAL RESIDUAL ERROR
!
      call ResidueCG (ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibp1, jbp1&
         , tiny, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, &
         c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, phi, &
         ex, ey, ez, vvol, dive, srce, rdt, scalsq, residu)
!
      rsqsum_new=0.0d0
      do n=1,ncells
         rsqsum_new = rsqsum_new+residu(ijkcell(n))*residu(ijkcell(n))
      enddo
      write(*,*)'PoissonCG:  dnorm after=',rsqsum_new
!
!
!
      return
      end subroutine Poisson

subroutine ResidueCG(ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibp1&
         , jbp1, tiny, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, &
         c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, &
         phi, gradphix, gradphiy, gradphiz, vvol, dive, srce, rdt, scalsq, &
         residu)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ncells
      integer  :: nvtx
      integer  :: iwid
      integer  :: jwid
      integer  :: kwid
      integer  :: ibp1
      integer  :: jbp1
      real(double)  :: tiny
      real(double) , intent(in) :: rdt
      real(double) , intent(in) :: scalsq
      integer  :: ijkcell(*)
      integer  :: ijkvtx(*)
      real(double)  :: c1x(*)
      real(double)  :: c2x(*)
      real(double)  :: c3x(*)
      real(double)  :: c4x(*)
      real(double)  :: c5x(*)
      real(double)  :: c6x(*)
      real(double)  :: c7x(*)
      real(double)  :: c8x(*)
      real(double)  :: c1y(*)
      real(double)  :: c2y(*)
      real(double)  :: c3y(*)
      real(double)  :: c4y(*)
      real(double)  :: c5y(*)
      real(double)  :: c6y(*)
      real(double)  :: c7y(*)
      real(double)  :: c8y(*)
      real(double)  :: c1z(*)
      real(double)  :: c2z(*)
      real(double)  :: c3z(*)
      real(double)  :: c4z(*)
      real(double)  :: c5z(*)
      real(double)  :: c6z(*)
      real(double)  :: c7z(*)
      real(double)  :: c8z(*)
      real(double)  :: vol(*)
      real(double)  :: phi(*)
      real(double)  :: gradphix(*)
      real(double)  :: gradphiy(*)
      real(double)  :: gradphiz(*)
      real(double) , intent(in) :: vvol(*)
      real(double)  :: dive(*)
      real(double) , intent(in) :: srce(*)
      real(double) , intent(out) :: residu(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
      real(double) :: rvvol
!-----------------------------------------------
!
!     calculate residue
!
      call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
         , c5z, c6z, c7z, c8z, phi, gradphix, gradphiy, gradphiz)
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, gradphix, gradphiy, gradphiz&
         )
!
      do n = 1, nvtx
         ijk = ijkvtx(n)
         rvvol = 1./vvol(ijk)
         gradphix(ijk) = -gradphix(ijk)*rvvol
         gradphiy(ijk) = -gradphiy(ijk)*rvvol
         gradphiz(ijk) = -gradphiz(ijk)*rvvol
      end do
!
!
!     calculate div(grad(phi))
!
      call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, vol, gradphix, gradphiy, gradphiz, dive)
!
!
!
      do n = 1, ncells
         ijk = ijkcell(n)
!
!
         residu(ijk) = (srce(ijk)-rdt*phi(ijk)-scalsq*dive(ijk))*vol(ijk)
!residu(ijk)=srce(ijk)-phi(ijk)
!
      end do
!
      return
      end subroutine ResidueCG

